Photonic band structures of periodic arrays of pores in a metallic host: tight-binding 

beyond the quasistatic approximation 



Kwangmoo Kim 1 ' 2 and D. Stroud 1 

1 Department of Physics, The Ohio State University, Columbus, Ohio 43210, USA 
2 School of Physics, Korea Institute for Advanced Study, Seoul, 130-722, Republic of Korea 

(Dated: November 11, 2011) 

We have calculated the photonic band structures of metallic inverse opals and of periodic linear 
chains of spherical pores in a metallic host, below a plasma frequency lj p . In both cases, we use a 
tight-binding approximation, assuming a Drude dielectric function for the metallic component, but 
without making the quasistatic approximation. The tight-binding modes are linear combinations of 
the single-cavity transverse magnetic (TM) modes. For the inverse-opal structures, the lowest modes 
are analogous to those constructed from the three degenerate atomic p-states in fee crystals. For the 
linear chains, in the limit of small spheres compared to a wavelength, the results are the "inverse" 
of the dispersion relation for metal spheres in an insulating host, as calculated by Brongersma et 
al. [Phys. Rev. B 62, R16356 (2000)]. Because the electromagnetic fields of these modes decay 
exponentially in the metal, there are no radiative losses, in contrast to the case of arrays of metallic 
spheres in air. We suggest that this tight-binding approach to photonic band structures of such 
metallic inverse materials may be a useful approach for studying photonic crystals containing metallic 
components, even beyond the quasistatic approximation. 



I. INTRODUCTION 

The photonic band structures of composite materials 
have been studied extensively. Such band structures are 
defined by the relation between frequency u> and Bloch 
vector k in media in which the dielectric constant is a 
periodic function of position. A major reason for such in- 
terest is the possibility of producing photonic band gaps, 
i.e., frequency regions, extending through all k-space, 
where electromagnetic waves cannot propagate through 
the medium. Such media have many potentially valuable 
applications, including possible use as filters and in films 
with rejection- wavelength tuning. 1] In systems with a 
complete photonic band gap, the spontaneous emission of 
atoms with level splitting within the gap can be strongly 
suppressed. @ 

Since light cannot travel through the photonic band 
gap materials (Bragg diffracted backwards), one of their 
applications can be a complete control over wasteful 
spontaneous emission in unwanted directions when a de- 
vice, such as a laser, is embedded inside a 3D photonic 
crystal. Q 2D photonic crystals can be used as optical mi- 
crocavities, microresonators,Q waveguides, [j| lasers, [f| 
or fibersQ while ID photonic crystals can be used as 
Bragg gratings or optical switches. [1] 

The photonic band structure of a range of materials 
has been studied using a plane wave expansion method. 
Typically, the method converges easily when the dielec- 
tric function is everywhere real, but more slowly, or not 
at all, when the dielectric function has a negative real 
part, as occurs when one component is metallic. For ex- 
ample, McGurn et a/.[9( used this method to calculate 
the photonic band structure of a square lattice of metal 
cylinders in two dimensions (2D) and of an fee lattice of 
metal spheres embedded in vacuum in 3D. They found 
that that method converged well when the filling fraction 
/ (i.e., volume fraction of metal spheres or cylinders) sat- 



isfied / < 0.1%. 

Kuzmiak et aZ.[l(j used the same method to calculate 
the photonic band structures for 2D metal cylinders in 
a square or triangular lattice in vacuum. For low / and 
cj > w p , the calculated photonic band structures are just 
slightly perturbed versions of the dispersion curves for 
electromagnetic waves in vacuum. However, for co < u p 
and H-polarized waves (magnetic field H parallel to the 
cylinders), they obtained many nearly flat bands for lj < 
w p ; these bands were found to converge very slowly with 
increasing numbers of plane waves. They later extended 
this work to systems with dissipation. ll| To describe 
dispersive and absorptive materials, they used a complex, 
position-dependent form of dielectric function. They also 
introduced a standard linearization technique to solve the 
resulting nonlinear eigenvalue problem. 

Zabel et aZ.[l2[ extended the plane wave method to 
treat periodic composites with anisotropic dielectric func- 
tions. In particular, they studied the photonic band 
structures of a periodic array of anisotropic dielectric 
spheres embedded in air. They found that the anisotropy 
split degenerate bands, and narrowed or even closed the 
band gaps. Much further work on anisotropic photonic 
materials has been carried out since this paper (sec, e.g., 
Ref. 0). 

A different type of periodic metal-insulator composite 
is a periodic arrangement of metallic spheres in an insu- 
lating host. Brongersma et aZ.[l3[ studied the dispersion 
relation for coupled plasmon modes in such a linear chain 
of equally spaced metal nanoparticles, using a near-field 
electromagnetic (EM) interaction between the particles 
in the dipole limit. They also studied the transport of 
EM energy around the corners and through tee junctions 
of the nanoparticle chain-array. 

Park and Stroud [3] also studied the surface-plasmon 
dispersion relations for a chain of metallic nanoparti- 
cles in an isotropic medium. They used a generalized 
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tight-binding calculations, including all multipoles. This 
approach is more exact than the previous point-dipole 
calculation, [l3j in a quasistatic limit, but still leaves out 
non-quasistatic effects associated with radiative damping 
(i.e., effects associated with the non- vanishing of V x E, 
where E is the electric field. They calculated the lowest 
bands as well as many higher bands and compared their 
results with those in Ref. 



13] 



Weber and Ford[l5| have shown that all calculations 
within the quasistatic approximation omit important in- 
teractions between transverse plasmon waves and free 
photon modes, even if the interparticle separation is small 
compared to the wavelength of light. Thus, most qua- 
sistatic calculations need to have certain corrections in- 
cluded at particular values of the wave vector. 

Recently, Gaillot et aZ.[ll| have studied the photonic 
band structures of another type of structure, a so-called 
inverse opal structure. This structure is an fee lattice of 
void spheres in a host of another material. Such a struc- 
ture can be prepared, e.g., starting from an opal structure 
made of spheres of a convenient substance, infiltrating it 
with another material, then dissolving away the spheres. 
In the work of Ref. [l6[ , the photonic band structure of Si 
inverse opal was calculated as a function of the infiltrated 
volume fraction / of air voids using three-dimensional fi- 
nite difference time domain (3D FDTD) method. It was 
found that for certain values of /, a complete band gap 
opens up between the eighth and ninth bands. 

In the present work, first we study the photonic band 
structure of an inverse opal structure, such as that in- 
vestigated in Ref. [l6[ , but instead of dielectric materials 
such as Si, we consider metals as the infiltrated materials. 
Thus, the material we study is also the inverse of the fee 
array of metal spheres studied by McGurn et al. [|[ Such 
metallic inverse opal structures have recently become of 
great interest, because it has been found that Pb inverse 
opals exhibit superconductivity. [l7[ These workers have 
studied the response of these materials to an applied mag- 
netic field, and have found a highly non-monotonic frac- 
tional flux penetration into the Pb spheres as a function 
of the applied field. 

As a second example, we study the photonic band 
structure of a linear chain of nanopores in a metallic 
medium. This is an inverse structure of a linear chain 
of metallic nanospheres, of which the dispersion relation 
is given in Ref. [l3j . As anticipated, we get a kind of 
"inverse image" of the dispersion relation found by Ref. 
13J in our system. 

For both types of structures, our primary method for 
studying the photonic band structures below the plasma 
frequency w p is a tight-binding approximation which 
is valid even in the non-quasistatic regime. Because 
the analogs of the tight-binding atomic states decay ex- 
ponentially in the metallic host medium, the resulting 
tight-binding waves do not lose energy radiatively, as do 
the corresponding waves along one-dimensional chains of 
metallic nanoparticles in air. Furthermore, because the 
modes are expanded in "atomic" states rather than plane 



waves, there is no convergence problem as there can be 
in the plane wave case. 

The remainder of this paper is organized as follows. 
In Section [IIJ we first present the formalism for calcu- 
lating the transverse magnetic (TM) and transverse elec- 
tric (TE) modes of a single spherical cavity in a metallic 
host. We then describe the method for calculating the 
photonic band structures of metallic inverse opals and of 
linear chains of nanopores in a metallic host, using a sim- 
ple tight-binding approach for lj < ui p . In Section UTT1 we 
give the numerical results for the TM and TE modes of 
a single cavity and those of the tight-binding method for 
the metal inverse opals and the linear chain of nanopores. 
Section [TV] presents a summary and discussion. 



II. FORMALISM 

In this section, we present a summary of the equations 
determining the band structure of a photonic crystal con- 
taining a metallic component with Drude dielectric func- 
tion e(ui) = 1 — oj'p/lo 2 and an insulating component of 
dielectric constant unity. The insulating component is 
assumed to be present in the form of identical spheri- 
cal cavities of radius R. We first write down the equa- 
tions for the TM and TE modes of a spherical cavity in a 
Drude metal. Then, we present a tight-binding method 
for uj < lu u . 



A. Spherical Cavity 

As a preliminary to calculating the photonic band 
structure, we first discuss the modes of a single spher- 
ical cavity in a Drude metal host. We begin with the 
TM modes of the cavity, then the TE modes. 



1. TM Modes 

It is convenient to describe the modes of the embedded 
cavity in terms of the B field. To that end, we combine 
the two homogeneous Maxwell equations 



V x E = — B, 



_ ,_. lid _ 

V x B = eE, 



to obtain a single equation for B 
1 „, 



V x 



1 - 6(x) > V x B 



1 



(1) 

(2) 
(3) 



Here, we have expressed the position- and frequency- 
dependent dielectric function e(x, uS) as l/e(x, u>) — 
0(x)/(l — Wp/w 2 ) + 1 — 0(x), where the step function 
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0(x) = 1 inside the metallic region and 0(x) = else- 
where. Multiplying this equation by cu 2 — lu 2 and simpli- 
fying, we obtain 

[u 2 - c 2 (l - 6>(x))]V x (V x B) = ^(c 2 - u, 2 )B. (4) 

Thus, inside the spherical void, we have 

, ,2 



V x (V x B) = — 



while inside the metal, 



V x (V x B) 



B. 



^B. 



(5) 



(6) 



For a spherical void within a metallic host, it is con- 
venient to solve these equations in spherical coordinates. 
It is readily found that, for the TM modes, the non- 
vanishing components of the solutions for B and E for 
Eq. © and Eq. © are [3 

B^r, 9) = ^lp} (cos 9), 



Ea.i-n — — 



-1(1- 



uir r 
ic dui(r 



l)-^P/(cosl 



P/(cos( 



cjr dr 

u e (r) = r[Aiji(kr) + B e n e (kr)], 



(7) 



where k — u/c, ji and ng are spherical Bessel functions, 
and the subscripts 0, r, and 9 denote components of the 
corresponding fields in spherical coordinates. 

Likewise, the solutions of Eqs. © and © within the 
metal will be 

= ^P/(cos0), 



En 



r 
ic 



ujr ( 1 



£(£ + l)^Wp £ ( COS ( 



ur ( 1 — 



dveir) P}(cose), 



dr 



ve{r) = r[C e j e (k'r) + D e n e {k'r)}, 



(8) 



where k' — J u 2 — lo 2 /c. 

The requirements that the normal displacement D and 
tangential E should be continuous at R gives the two 
conditions 



and 

due(r) 



di 



u t {R) = v t (R) 
dvi(r) 



r=R 



1 — -E- 



di 



k 2 dvi{r) 



r=R 



k' 2 dr 



(9) 



r=R 

(10) 



where R is the radius of the spherical cavity. Since the 
fields at the center of the void sphere must be finite, we 
also have 



ui(r) = rjt{kr), 



(11) 



where we have normalized the solution so that the coef- 
ficient At = 1. From Eq. © we have 

MkR) = C e je{k'R) + D e n e (k'R), (12) 

whereas, from Eq. (ITU)) , we get 

k' 2 \jt{kR) + kRj' f (kR)] 
= k 2 [C t j t {k'R) + D e n e (k'R) + k'R{Cd' e (k'R) 
+ D e n' e {k'R)}} . (13) 

The coefficients Ci and Di can then be determined from 
these boundary conditions. 

The results for uj < w p , can be obtained by making 
the substitution k' — > ik' , with k' real. In this case, the 
radial component of Eq. © takes the form 



d 2 u e 
dr 2 



2^P 



ut = 0, 



(14) 



which is the modified Bessel equation. The solutions of 
Eq. (TT4")) are the modified spherical Bessel functions ie and 
kg (note that this kg is different from the wave vectors k 
and k'). In this case, Eq. © becomes 



v e (r) = r[C e i e {k'r) + D e k e (k'r)}, 



where k' = ^Jlu 2 — lu 2 /c. In addition, Eqs. ([10 
and p3p are transformed, respectively, into 



(15) 



du e {r) 



dr 



r=R. 



k 2 dvi(r) 
k' 2 dr 



r=R 



MkR) = Cm{k'R) + D e k e {k'R), 



(16) 



(17) 



and 



k' 2 [j e (kR) + kRj' e (kR)] = -k 2 [CMk'R) 
+D e k e {k'R) + k'R{Cii' e {k'R) + D e k' e (k'R)}] .(18) 

It is of interest to consider the specific case of a spher- 
ical cavity in an infinite medium. In this case, Ci = 
because ii(x) diverges at large x. As a result, Eqs. (|15|) . 
(Tl~7| . and (fTS]) become, respectively, 



vg(r) = rDiki{k r), 



j e (kR) = D e k e (k'R), 



and 



k' 2 [MkR) + kRj' £ (kR)} = -k 2 D e [k e {k'R) 
+k'Rk' ( ,(k / R)} . 



(19) 



(20) 



(21) 
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From Eq. (|20| we have Di — je(kR)/ke(k' R), and hence 
Eq. (j2"Tj) becomes 



fc' 2 {j e (kR) + kRji(kR)} = -fc 2 ^|^ Mk'R) 



+k' Rk' e {k' R)] . 



(22) 



We can readily obtain the asymptotic forms of the solu- 
tions when kR <C 1 and k'R <C 1. In this case ji{kR) and 
ke(k'R) have the asymptotic forms ji(kR) m (kR) e / (2£ + 
1)!! and k e (k'R) w v^T(£ + l/2)2*~ V(fc'i?f +1 . In this 
limit, Eq. (f2"2"j) . after some algebra, reduces to simply 



k' 2 (£ + 1) = fc 2 £. 



Since k' = Ji>j 2 — w 2 /c, Eq. (|23|) is equivalent to 



2 1 + 1 2 

21+1 p 



(23) 



(24) 



The largest value, u; = y/2/ 3w p , occurs at £ = 1 and the 
limiting value for large £ is to — <jj p /y2. 



2. TE Modes 



For the TE mode, inside the spherical void, we have 



V x (V x E) = — -E, 
whereas inside the metal, we have 



V x (V x E) = 



(25) 



(26) 



We now use these equations to calculate E^, B r , and 
Bg. From Eq. (|2"5j) we get 

i^r,0) = — ^ (cos 0), 



B 



ue(r) 



,,„ — £(£+l)^^P,(cos0), 



ic <9i^(r) . 



P/(cos( 



in — .. 

wr or 

u/(r) = r[A^(fcr) + B e n e (kr)}, (27) 
where k = lu/c. The solutions of Eq. (|26| arc 

^,out(r,0) = Pi {cos 9), 

B r ,out = — e(t + i)^P t (poa6), 
tor r 

Be, out = jr+Pl (cos 9 ), 

lot or 

ve(r) = r[C e ji(k'r) + D m (k'r)}, (28) 



where 



k! = ^Jijj 2 — ujp/c. 



Since normal B and tangential H should be continuous 
on the boundaries, we obtain the conditions 



ui(R) = vt(R), 
as in the TM case, and 



dui(r) 



Or 



dvg(r) 



r=R 



Or 



(29) 



(30) 



r=R 



Since the fields must be finite at the center of the void 
sphere, we can choose 



u e (r) = rji(kr), 



(31) 



we also take the coefficient Ag — 1. From Eq. (|29p we 
have 



U (kR) = Gatik'R) + D e n e (k'R), 
while, from Eq. (j3"0")l . we get 



(32) 



u(kR) + kRj[{kR) 
= C t u(k'R) + D e n e (k'R) + k'R [CijfflR) 

+D e n' e {k'R)}. (33) 

The corresponding equation for u> < u> p , can again be 
obtained by the transformation k' — > ik! . The TE modes 
for w < w p using the modified spherical Bessel functions 
it(x) and k e (x) satisfy Eqs. ((Til), (15]), ([30]), and (fT7| . 
The only changes are in Eq. ([T8|) . which becomes 



j e (kR) + kRj' t {kR) = C e i e (k'R) + D e k e (k'R) 
+k'R{Cii' e (k'R) + D e k' e (k'R)}. (34) 

In an infinite medium, these conditions become, from Eq. 

u(kR) + kRj' t {kR) = /4!S Mk'R) + k'Rk' e (k'R)} . 

If we consider the asymptotic forms of the solutions 
when kR -C 1 and k'R < 1 as we did for the TM modes, 
Eq. ([35]) simplifies to 



1 



(36) 



which gives £ = —1/2. Since £ must be a positive integer, 
we see that there are no eigenvalues for TE modes in the 
limit kR < 1 and k'R < 1. 



B. Tight-Binding Approach to Modes for u> < oj p 

We now turn from describing the single-cavity modes 
to a discussion of the band structure for a periodic array 
of such cavities. In conventional periodic solids, the tight- 
binding method is very useful in treating narrow bands. 
In what follows, we try to suggest an analogous tight- 
binding approach for the lowest set of TM modes in a 



5 



periodic lattice of spherical cavities in a metallic host, 
in the frequency range cj < ui p . We apply the resulting 
method, first, to an fee lattice of pores, and then to a 
linear chain of spherical pores in a metallic host. 

Even though these are TM modes, it is convenient to 
describe them now in terms of their electric fields. We 
denote the electric field of the Ath mode by Ea (x) . This 
field satisfies 



Vx(VxE A (x))- 



, E A (x) = OEa(x) = -£E A (x), 

(37) 

where O = V x (Vx) + (u; 2 /c 2 )6>(x) is the "Hamilto- 
nian" of this system. Since O is a Hermitian opera- 
tor, the eigenstates corresponding to unequal eigenvalues 
w 2 /c 2 and w 2 /c 2 are orthogonal and may be chosen to 
be orthonormal. (The orthogonality may also be proved 
directly by integration by parts.) The orthonormality 
relation is 



E A (x) • E„(x)dx = S 



(38) 



Since Ea(x) is real for u> < cj p , the complex conjugation 
is, in fact, unnecessary. 

In Sec. Ill A 11 our paper already gives the equations 
determining the electric and magnetic fields of isolated 
TM modes for a spherical cavity. The lowest set corre- 
sponds to i = 1, and there should be three of these. For 
a spherical cavity, all three are degenerate, i.e., all three 
have the same eigenfrcquencies. Even though the three 
modes have equal frequencies, one can always choose an 
orthonormal set, with electric fields Ei, E2, and E3 sat- 
isfying the orthonormality relation in Eq. (|38l) . 

In order to obtain the tight-binding band structure 
built from these three modes, we need to calculate matrix 
elements of the form 



E*(x) -OE^x-R^x, 



(39) 



corresponding to two single-cavity modes associated with 
different cavities centered at the origin and at R. Here, O 
is the "Hamiltonian" of the system as defined implicitly 
in Eq. (|57)l . 

Next, we introduce normalized Bloch states associated 
with the three 1=1 single-cavity modes. In order to 
do this, we first make the usual tight-binding assump- 
tion that the "atomic" states corresponding to different 
cavities are orthogonal: 



J E A (x - R) ■ E p (x - R')dx - 8 X ,„S R>RI 



(40) 



This orthogonality of states on different cavities is reason- 
able since the fields fall off exponentially with separation. 
The orthonormal Bloch states then take the form 

E k , A (x) = N- 1 ' 2 J2 e ik - R E A (x - R), (41) 

R 



where k is a Bloch vector, and the R's are the Bravais 
lattice vectors. In writing Eq. ([4"Tj). we have assumed 
that there are N identical spherical cavities, and that 
the Bloch states satisfy the usual periodic boundary con- 
ditions of Born- von Karman type. We also introduce the 
elements of the "Hamiltonian" matrix 



M A , M (k) =^e lk - R AfA iAi (R). 



(42) 



R 



We can then obtain the frequencies w(k) by diagonal- 
izing a 3 x 3 matrix as follows: 



det 



( k ) gat 1 s 

c 2 c 2 



A . / / 



= 0, 



(43) 



where w a t is the eigenvalue of a single-cavity mode. The 
solutions to these equations give the three p-bands for 
a periodic lattice of cavities in a metallic host. This 
procedure is analogous to that used in the well-known 
procedure for obtaining tight-binding bands from three 
degenerate p-bands in the electronic structure of conven- 
tional solids (see, for example, Ref. [19]). 

We briefly comment on the connection between this ap- 
proach and that used by earlier workers. [HI [l4| In this 
work, the authors treat wave propagation along a chain 
of metallic nanoparticles. They use the tight-binding ap- 
proximation, as we do, but in the quasistatic approxi- 
mation in which one assumes that V x E = 0. This 
approximation is reasonable when both the particle radii 
and the interparticle separations are small compared to 
a wavelength, but is not accurate in other circumstances. 
Furthermore, even in the small-particle and small-sepa- 
ration regime, this approximation still fails to account 
for the radiation which occurs at certain wave numbers 
and frequencies. The present approach would generalize 
this tight-binding method to (a) three dimensions as well 
as one; (b) pore modes instead of small particle modes; 
and most importantly (c) larger pores and larger inter- 
particle separations, via extension beyond the quasistatic 
approximation. 

Next, we discuss the numerical evaluation of the re- 
quired matrix elements, Eq. ((391) . The relevant electric 
fields are given in this paper, but in spherical coordinates. 
It is not difficult to convert these into Cartesian coordi- 
nates. The operator O is just a little trickier. We first 
note that O — Or + O', where Or is the single-cavity 
operator: Or = V x (Vx) if x is inside the Rth cavity 
and Or = V x (Vx) + uj 2 /c 2 otherwise. Now we also 
have 



OrE^x-R) = ^E /3 (x- 



(44) 



since is an eigenstate of Or with an eigenvalue w 2 t /c 2 . 

But since we are assuming that the overlap integral be- 
tween "atomic" electric field states centered on different 
sites vanishes, the term involving Or does not contribute 
to the matrix element M a ^, which is therefore just given 
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by 



M aJj (R) = J E Q (x) • C'E^(x - R)dx. (45) 



We can also write 



R' 



where 



?R' 



(x-R'), 



(46) 



(47) 



is a step function which is unity inside the cavity centered 
at R' and is zero otherwise. 

A reasonable approximation to Eq. (|46l) might be to 
include just R' = 0. In this case, we finally will get 

M a>0 (R) ~ ^ J E a (x) • E^(x - R)dx, (48) 

where the integral runs just over the cavity centered at 
the origin. As a further approximation, we can just re- 
place E / g(x— R) by the value of this function at the origin, 
i.e., Ep(— R). Then this field can be taken outside the 
integral and we just have 

M Qi/? (R) ~ ^E/j(-R) ■ J E Q (x)dx, (49) 

where once again the integral runs over the cavity cen- 
tered at the origin. 

Next, we attempt to calculate the relevant quantities 
needed to solve for this matrix element. In order to use 
the tight-binding approach we will need to normalize the 
individual eigenstates E a . Therefore, we will begin by 
obtaining this normalization. For I = 1, the M^(r)'s are 
r times spherical Bessel functions. We write this field as 



E, 



2d 
kr 



-r-, — ji(fcr)cos( 



(50) 



where we have used the relation Pi (cos 9) = cos 9, and 
introduced the normalization constant C\, which will be 
determined below. Similarly, 



En 



C x d[rh{kr)] 



sin 9, 



(51) 



' ln kr dr 
where we use Pl(cos9) — — sin 9. For r > R, we have 
2D C 

£r,out = +—p-±k 1 (tfr)coB6 (52) 

and 

n.Hi r)\rh.^ (h'rW 

(53) 



Dtddyrktik'r)] 
Ee.out = sm f 



k'r dr 

We will need the integrals of the Cartesian components 
of the field over the volume of the sphere centered at the 



origin. Let us assume we are considering the z mode, 
i.e., the one for which 9 refers to the angle from the z 
axis. Then the symmetry of the problem shows that only 
the z component of the electric field will have a nonzero 
integral. Also, we have that 



E z ,i n (r, 9) = E r - m cos 9 - Eg >iD sin ( 



(54) 



Thus, after a little algebra, we find that the integral of 
this field over the volume of the cavity is 



E z ,- m (r, 9)dv 



8ttCi 
3fc 



R 2 3i(kR). 



(55) 



Next, we work out the coefficient D\. It is determined 
by the boundary conditions at r — R. These conditions 
are that D r and Eg should be continuous at r = R. These 
two conditions determine not only the value of D\ but 
also the allowed frequency. After a bit of algebra, we find 
that 



D, 



uj 2 jx{kR) k nikR) 



c 2 kk' ki(k'R) k'k x {k'R)' 



(56) 



The allowed value of lu is given by Eq. (f2"2"j). 

Finally, we need the normalization constant C\. We 
choose this so that the integral of the square of the elec- 
tric field for a single-cavity mode should be normalized 
to unity. This condition may be written 



2C\ 4tt 
2C\D\ 4tt 



2jt(kr) 



^ d[rjx(kr)] 



dr 



2k 2 (k'r) 



d[rki(k'r)] 
dr 



= 1. 
If we write 



dr 



(57) 







and 



k'R 



2j 2 (x) + 



2k 2 1 (x) + 



d[xji(x)] 
dx 



d[xki(x)] 
dx 



dx^FxikR) (58) 



dx = F 2 {k'R) 1 (59) 



then we can express the normalization condition as 



8irC 2 



F^kR) 
k 3 



Di 



F 2 {k'R) 



k' 3 



(60) 



Therefore, we can now write out an explicit expression 
for the matrix element M Qj/ 3(R) given in Eq. (|4l?|). For 
the ath mode, the integral of E Q over the volume of a 
cavity is a vector in the ath direction. To evaluate Eq. 
(14T)1) , we need the component of the ath mode in the /3th 
direction at a position R. Let us first consider the zth 
mode (a — z). We can use Eqs. (fS"2")) and (|5"3")l to rewrite 



this field in Cartesian coordinates with the additional 
equations 



cos 9 



sinf 



\/x 2 + y 2 



xx + yy + zz 



xzx + yzy \/ x 2 + ; 
ry/ x 2 + y 2 r 



(61) 



We just substitute these expressions back into Eqs. (j5"2")l 
and (155)) to get the Cartesian components of the field for 
a mode parallel to the z axis. For the mode parallel to 
the x axis, we just permute the coordinates cyclically: 
z — > x, x — > y, and y z. Similarly, for the y modes, we 
make the permutation {x,y,z) — > {z,x,y). 

Using these results, we should be able to compute all 
the elements in the tight-binding matrix and hence obtain 
the band structure for the photonic p-bands in the tight- 
binding approximation, in either one or three dimensions. 



III. NUMERICAL RESULTS 

For the inverse opals we arbitrarily assume a lattice 
constant d = 500-\/2 nm, and a void sphere radius R = 
150 nm as in Fig. HJa). This choice is the same as that 
of Ref. [17j . where the Pb inverse opal has this lattice 
constant. Since the volume of the primitive unit cell is 
v c = d 3 /4, this corresponds to a void volume fraction 
/ = 0.160. For the linear chain of nanopores (see below) 
this d is the separation between two nanopores and R is 
the radius of a nanopore as in Fig. [TJb). 

Our band structures for the inverse opals are expressed 
in terms of the standard notation for k values at sym- 
metry points in the Brillouin zone. These are T = 
(0,0,0), X = (27r/<2) (0,0,1), U = (27r/d) (1/4, 1/4,1), 
L = (27r/d)(l/2,l/2,l/2), W = (2vr/d)(l/2, 0, 1), and 
K = (27r/d)(3/4, 0,3/4). 

The metallic dielectric functions we assume for the in- 
verse opals and linear chain of nanopores are of the usual 
Drude form, 




FIG. 1: (Color online) Schematic diagram for (a) an inverse 
opal structure with a lattice constant d and a void sphere ra- 
dius 7?; (b) a linear chain of nanopores with a pore separation 
d and a nanopore radius R. 



wave vector kd, and the band structures are parameter- 
ized by the two constants oj p d/c and / for the case of 
inverse opals. 



e(w) = 1 



_P 

o ) 



(62) 



where uj p is the plasma frequency of the conduction elec- 
trons. e(uj) < when uj < uj pi while e(w) > when 
uj > uj p . Our calculations are thus carried out assuming 
that the Drude relaxation time r — > oo. For a metal in its 
normal state, lo 2 = 4iTne 2 /m, where n is the conduction 
electron density and m is the electron mass. Note that 
with this choice of dielectric function, the entire band 
structure can be expressed in scaled form. That is, the 
scaled frequency ujdjc is a function only of the scaled 



Since we are considering void spheres in inverse opals 
and linear chains of nanopores, it is of interest to consider 
electromagnetic wave modes in a single cavity, which 
could be considered a single "atom" of the void lat- 
tice. We show only results for oj < uj p , since these are 
the results most relevant to possible narrow-band pho- 
tonic states in the inverse opal structure. Our results 
for oj < UJ P for an isolated spherical cavity in an infinite 
medium, and when kR <s^. 1 and k' R <C 1 are given in 
Table |U These two inequalities are reasonable for our in- 
verse opal system parameters d — 500\/2nm, R = 150nm, 
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and uJ v d/c = 1, because 



0.2121, 



k'R = 




= v/0.045 - (kR) 2 < V0.045 = 0.2121 



The (modified) spherical Bessel functions in Eq. (|2"2")l arc 
extremely close to the ui axis for £ > 5, so that it is 
difficult to get eigenfrequencies for I > 5 in the isolated 
spherical cavity. However the eigenfrequencies continue 
to exist even for I > 5 when kR <C 1 and k'R <C 1. 





Infinite medium 


kR < 1, k'R « 1 


1 = 


1 


0.1296 


0.1299 


£ = 


2 


0.1232 


0.1233 


1 = 


3 


0.1203 


0.1203 


£ = 


4 


0.1186 


0.1186 


1 = 


5 


0.1178 


0.1175 



TABLE I: TM mode frequencies to' — Lud/(2nc), where lj < 
LOp and ujpd/c — 1, calculated for an isolated spherical cavity 
("Infinite medium") and those when both kR <C 1 and k' R -C 
1. The (modified) spherical Bessel functions are extremely 
close to the u/ axis for £ > 5, so that it is difficult to get 
eigenfrequencies for £ > 5 in the isolated spherical cavity. 
However this does not happen when kR <C 1 and k'R <C 1. 



The solutions to Eq. (1551) do not exist for uj < uj p 
with ujpd/c = 1. This fact is consistent with that the 
eigenvalues for uj < uj p do not exist for TE modes when 
kR < 1 and k'R < 1. 

For our fee calculations, we calculate the band 
structure including only the 12 nearest-neighbors of 
the cavity at the origin. Thus R = (d/2)(±l, ±1, 0), 
(d/2)(±l, T l,0), (d/2)(±l,0,±l), (d/2)(±l,0,Tl), 
(d/2)(0,±l,±l), and (d/2)(0,±l,Tl). Assuming 
Wpd/c = 1 and using w at <i/(27rc) = 0.1296 for £ = 1 
in an infinite medium, we get the tight-binding results 
in Fig. [3J This figure shows three separate bands in 
the X-U-L region and X-W-K region as expected 
for the p-bands. The bandwidth is relatively small as 
M ajj a(R)eP ~ 0.001, which proves the general relation 
between the bandwidth and the overlap integral. fl9j 
All three bands are degenerate at k = (the T point). 
In addition, there is a double degeneracy when k is 
directed along either a cube axis (T-X) or a cube body 
diagonal (T-L), the higher (concave upward) bands 
being degenerate in both cases. The lower two bands 
have a band gap at the U point, and these bands cross 
at the W point. 

Next, we turn to the band structure of a periodic 
linear chain of spherical nanopores in a Drude metal 
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FIG. 2: Tight-binding inverse opal band structure for uj < uj p 
with d = 500\/2 nm, R — 150 nm, and uj p d/c = 1.0, using 
ujaxd/ {2nc) — 0.1296 for £ — 1 in an infinite medium. The 
horizontal dotted line represents the "atomic" level. 



host. For this linear chain, the Bravais lattice vectors are 
R = d(0, 0,±n), where ±n is the nth nearest-neighbor, 
d is the separation between two nanopores and we as- 
sume that the chain is directed along the z axis. We can 
calculate the tight-binding band structure including as 
many sets of neighbors ±n as we wish. To compare our 
results with those in Ref. [l3j], we first use their param- 
eters, R — 25 nm and d = 75 nm, together with their 
overlap parameter u)\ — 1.4 x 10 15 rad/s. These com- 
bine to give io-pd/c = 0.35. The "atomic" frequency is 
found by solving Eq. ([2"2"| and gives uj at d/ (2irc) = 0.0454 
for I = 1 in an infinite medium. Our resulting tight- 
binding dispersion relations are shown in Fig. [3] with 
only nearest-neighbors included. Note that our frequen- 
cies are given in unit of 2nc/d while the results of Ref. 
(l3| are not scaled. Our results are exactly the inverse 
images of theirs — that is, we would get their curves 
(to within a constant of proportionality) if we reflect our 
curves through the horizontal line of the atomic level, 
and the transverse (T) branches are twofold degenerate 
as are theirs, while the longitudinal (L) branch is non- 
degenerate. Our eigenfrequency for a single-cavity cj at 
corresponds to their resonance frequency ujq. As we in- 
crease the number of nearest-neighbors (nn's) included, 
the separation between the L and T branches increases 
at the zone center but decreases at the zone boundary, 
as shown in Fig. the same trend is seen in Fig. 1 in 
Ref. 13(. The sum also converges quickly, so there is 



little difference between the dispersion relation including 
through the next-nearest-neighbors and that including 
through the 5th nearest-neighbors. 

We have carried out similar calculations using other 
values of the parameter u p d/c 1 namely 1.0, 2.0, and 5.0. 
Such calculations are possible here because our calcu- 
lations are non-quasistatic, so that the overlap integral 
between neighboring spheres falls off exponentially with 
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FIG. 3: Tight-binding results of a periodic chain of nanopores 
in a Drude metal host, for uj < cj p . We use d = 75 nm, 
R — 25 nm, and u) p d/c = 0.35, using uj a td/(2nc) = 0.0454 for 
£ — 1 in an infinite medium. Only the nearest-neighbors are 
included. The horizontal dotted line represents the "atomic" 
level. In this and the following all plots, "L" and "T" denote 
the longitudinal and transverse branches, respectively. 
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FIG. 4: (Color online) Same as Fig. [3] but including three 
different numbers of neighbors: nearest-neighbors (nn's), 
next-nearest-neighbors (nnn's), and fifth-nearest-neighbors 
(5nn's). 



separation. The results are given in Figs. El and 
respectively. The corresponding results including more 
overlap integrals are shown in Figs. andfTUl respec- 
tively. It is also striking that, as uj p d/c increases in going 
from Fig.0to Figs. and0 the ratio r LT of the width 
of the L band to that of the T band steadily decreases. 
In Fig.0 r LT > 1, in Fig.0 r LT < 1, while in Fig. [7] (for 
which uipd/c — 2.0), 7*lt ~ 1- 

One could also say that, except for an overall scale fac- 
tor, Fig. [9] looks like an inverted image of Fig. about 
the horizontal line of w a t. The dispersion relations for the 




intermediate value oj p d/c = 2.0 has nearly perfect sym- 
metry about the horizontal line of w at (the T branches 
are nearly reflections of the L branch about the horizon- 
tal line of w at ), as in Fig. [7J For the nn case, the T 
and L bands cross at ±7r/(2d), as can be seen in Figs. 
El and When further neighbors are included, they 
cross at smaller values than |7r/(2d)|, as can be seen in 
Figs. and ITUI but the crossing points get closer to 
±7r/(2d) as w p increases. Also, the effects of including 
further neighbors become smaller as lu p increases; they 
are smallest at ui p d/c = 5.0, as can be seen in Fig. [TUl 




nfd 



FIG. 5: Same as Fig. except iv p d/c — 1.0 and w at d/(27rc) = 
0.1291. 




FIG. 6: (Color online) Same as Fig. [H except u} p d/c 
and uj at d/{2%c) = 0.1291. 



1.0 



Next we consider values of R/d other than 1/3, but 
still keeping the same value of wi = 1.4 x 10 15 rad/s (i.e., 
Wpd/c = 0.35). For a smaller R/d = 0.25, the varia- 
tion of the band energies with k becomes smaller, as seen 
in Fig. [TTJ than it is in Fig. but the crossing points 
between the L and T branches still occur at ±w/(2d). 



10 





0.262 r 




U.ZOU 




U.Zjo 


o 


n 

U.ZJO 


(N 




-a 


U.ZJ4 


3 


0.252 - 




0.250 - 




0.248 - 




0.246 ^ 






0.580 r 




0.578 ' : 




0.576 - 




0.574 - 


g 0.572 
5 0.570 
§0.568 
¥ 0.566 




0.564 ,. 




0.562 




0.560 




0.558 - 




-Ttid 



FIG. 7: Same as Fig. |3j except uu p d/c — 2.0 and u} a _ t d/(2nc) = 
0.2537. 



FIG. 9: Same as Fig. O except ui p d/c — 5.0 and w a td/ '(2nc) 
0.5691. 
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FIG. 8: (Color online) Same as Fig. |4j except tj p d/c = 2.0 
and u) a , t d/{2-Kc) = 0.2537. 



FIG. 10: (Color online) Same as Fig. 21 except u) p d/c 
and Lu^d/{2-irc) = 0.5691. 



5.0 



This behavior becomes clearer when the results for sev- 
eral values of R/d are plotted together as in Fig. [12J As 
R/d increases, the variation of the band energies with k, 
and the separation between the L and T branches at both 
the zone center and zone boundary, increase, but the L 
and T branches still cross at ±7r/ (2d). If we include more 
neighbors up to fifth nearest-neighbors, but consider only 
up to R/d — 0.4, we get the dispersion relations shown 
in Fig. Q2J These show the same trends as in Fig. [IJJJ ex- 
cept that the band crossing points occur at values of |fc| 
slightly less than |7r/(2d)|. Furthermore, the separation 
between the L and T bands increases slightly at k = 0, 
but decreases slightly at k = ±n/d. We show only R/d 
up to 0.4 in this Figure because, in the quasistatic limit, 
there is evidence that for larger values of R/d the disper- 
sion relations are significantly modified by higher values 
oft Gil 



IV. DISCUSSION 

In this work we have calculated the photonic band 
structures of metal inverse opals and of a linear chain 
of spherical voids in a metallic host for frequencies be- 
low w p , when t = 1 using a tight-binding approxima- 
tion. In both cases, we include only the £ = 1 "atomic" 
states of the voids. As a possible point of comparison, we 
have also computed the same band structures using the 
asymptotic forms of the spherical and modified spherical 
Bessel functions for small void radius. In this asymptotic 
region, there are only TM modes. The results for the lin- 
ear chain of voids can be considered as the "inversions" of 
those in Ref. [l3j], in the sense discussed earlier. In other 
words, if we reflect our L and T branches with respect 
to the atomic energy level, we would get their L and T 
modes. 
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FIG. 11: Same as Fig. [3] except that R/d = 0.25 and 
u;atd/(27rc) = 0.04546. 
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FIG. 12: (Color online) Plotting together six different re- 
sults for uj < Up, all with u p d/c = 0.35, but with different 
R/d: R/d = 0.33 and w at d/(27rc) = 0.04544; R/d = 0.40 
and cj at d/(27rc) = 0.045426; R/d = 0.45 and u at d/(2nc) = 
0.045412; fl/d = 0.49 and u at d/(27rc) = 0.045399; and 
R/d = 0.50 and a; at d/(27rc) = 0.045396. In each case, only 
nearest-neighbor overlaps are included. 



Although we did not discuss this approach, we did at- 
tempt to use the plane wave expansion method to cal- 
culate the band structure for the inverse opals, similarly 
to Refs. Q and [l(|. Just as found in those papers, the 
photonic bands for modes below w p depend on the num- 
ber of plane waves included in the expansion and on the 
type of field, B or E, used in the expansion. Further- 
more, this plane wave expansion method gives a large 
number of flat bands below u; p , which arc difficult to in- 
terpret physically. Because of this problem, and because 
of the apparent non-convergence of this approach with 
the number of plane waves, we do not present these re- 
sults here. By contrast, the band structures above ui v , 



FIG. 13: (Color online) Same as Fig. 1121 except that only 
three R/d's are plotted, with inclusion of up to the fifth 
nearest-neighbors. We omit the three largest values of R/d 
because it may be necessary to include more than just 1=1 
when R/d > 0.4. 



when calculated using the plane wave expansion, varied 
smoothly with k and converged well with the number of 
plane waves included. 

In calculating the tight-binding band structure for the 
linear chain (and the inverse opal structure), one should, 
in principle, include all the neighbors. But in practice, 
for the linear chain, it is sufficient to include only up 
to the fifth nearest-neighbors. This calculation is easily 
carried out, since the matrix in Eq. is already diag- 
onal in x, y, z and the sum converges quickly. In fact, 
even the inclusion of neighbors beyond the first two sets 
changes the band structure very little. The smallness of 
the further neighbor effect is particularly apparent when 
ujpd/c = 5.0, as in Fig. [101 

Although we studied the 3D and ID lattices, we have 
not investigated a 2D lattice of spherical pores in a metal- 
lic host. For the 2D case, the matrix element M ai p(R) 
can again be readily calculated using our tight-binding 
approximation. It is expected to have some nonzero off- 
diagonal elements in addition to the diagonal elements. 
Thus band structures somewhat similar to the 3D band 
structures shown in Fig. [2] are also expected in the 2D 
case. 

In the quasistatic case, for metal grains in air, when 
R/d is greater than about 0.4, it becomes important to 
include more than just I = 1 as in Ref. Inclusion 
of such higher £'s might be rather difficult in the present 
dynamical case, though it would be possible in the qua- 
sistatic limit for ID chains of spherical nanopores. 

In the present work, we have considered only the case 
of one cavity per primitive cell and t = 1 . It would be of 
a great interest to consider multiple cavities per unit cell. 
Of course, in this case, the dimension of the matrix in Eq. 
(f39|) will increase and the Bloch states in Eq. (|4"Tj) will 
acquire an additional index. It should be straightforward 
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to extend the present work to such a case, which would 
make an interesting subject for future work. 

The surface effect of the metal nanopores on the eigen- 
states becomes prominent as the radius of a pore (i?) or 
the ratio of radius to centcr-to-ccntcr separation (R/d) 
increases. Shockley surface states form on metal sur- 
faces depending on the type of metal, such as the Fermi 
wavevector fcp or the Fermi energy cf- These surface 
states will interact with the eigenstates, resulting in the 
change of eigenvalues. However we think these effects are 
negligible, especially in the quasistatic limit since in the 
asymptotic limit, these oscillatory interactions take the 
form|20| 

pv-ymr.s „ f 4e F \ sin(2fc F a + 29) 

where a is the interatomic separation, is the effective 
interaction phase shift, fcp the Fermi wavevector of the 
isotropic surface state, and e F the Fermi energy. The 
proportionality constant gives the consequences of scat- 
tering into bulk states. pO] The very slow a -2 decay of 
these interactions allows them to play a role at large sep- 
arations, but the overall magnitude is small due to the 
sinusoidal term. 

In summary, we have described a tight-binding method 



for calculating the photonic band structure of a periodic 
composite of spherical pores in a metallic host, and have 
applied it to both ID and 3D systems. The method is 
fully dynamical, and is not limited to very small pores. 
The method does not have the convergence problems 
found when the magnetic or electric field is expanded in 
plane waves. Furthermore, there are no radiation losses 
to consider, unlike the complementary case of small metal 
particles in an insulating host, because the fields associ- 
ated with these modes outside the pores are exponentially 
decaying. Thus, this method may be useful for a variety 
of periodic metal-insulator composites. It would be of 
interest to compare these calculations to experiments on 
such materials. 
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